% Determine initial pixel states
function [Lambda] = lambda(X)

% Init label vector based on pixel gray level
threshold = max(max(X))/2;
Lambda = uint8(X > threshold);

% Combine pixel states with neighbours' states
[M, N] = size(X);
for m = 2 : M-1
    for n = 2 : N-1
        % q*2^0 + r*2^1 + s*2^2 + t*2^3 + u*2^4
        Lambda(m, n) = bitand(Lambda(m, n), 1) ...
            + bitand(Lambda(m, n-1), 1) .* 2 ...
            + bitand(Lambda(m-1, n), 1) .* 4; ...
            + bitand(Lambda(m, n+1), 1) .* 8 ...
            + bitand(Lambda(m+1, n), 1) .* 16;
    end
end